home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / PARSER / KPARS_00 / KMATH.PAS < prev    next >
Pascal/Delphi Source File  |  1993-09-02  |  15KB  |  564 lines
  1. {$A-}
  2. {$B-}
  3. {$D-}
  4. {$E-}
  5. {$F+}
  6. {$I-}
  7. {$L-}
  8. {$N-}
  9. {$O+}
  10. {$R-}
  11. {$S-}
  12. {$V-}
  13.  
  14. UNIT KMath;
  15. {+H
  16. ---------------------------------------------------------------------------
  17.   Version     - 0.21
  18.  
  19.   File        - KMath.PAS
  20.  
  21.   Copyright (c) Klingon Software Services 1987..1993 except where noted.
  22.                 All rights reserved.
  23.  
  24.   Author      - Keith S. Brown (except where otherwise noted)
  25.                 Surface mail:              Email:
  26.                   K.Brown
  27.                   2437 Bay Area Blvd #20
  28.                   Houston, TX 77058 (USA)  Voice:713-486-6765
  29.  
  30.   Purpose     - Standard Mathematics Functions and Procedures.
  31.  
  32.   Remarks     - None.
  33.  
  34.   Language    - Borland International's Turbo Pascal V:4.x+ for MS-DOS
  35.  
  36.   Written     - 1988.0202 V:0.00 (KSB) Wrote initial version.
  37.   Revised     - 1988.1026 V:0.10 (KSB) Fixed ArcCos.
  38.               - 1990.0831 V:0.20 (KSB) Added Approx function
  39.               - 1990.0906 V:0.21 (KSB) Added Normalize procedure.
  40. ---------------------------------------------------------------------------}
  41. INTERFACE
  42.  
  43.  
  44. {}FUNCTION  Approx(a,b:REAL):REAL;
  45.     { Rounds b if Frac(b) > a; else truncs b
  46.     }
  47. {}FUNCTION  ArcCos(x:REAL):REAL;
  48.     { Returns arcCosine of a value between -1 and 1
  49.     }
  50. {}FUNCTION  ArcSin(x:REAL):REAL;
  51.     { Returns arcSine of a value between -1 and 1
  52.     }
  53. {}FUNCTION  ArcTan2(y,x:REAL):REAL;
  54.     { Returns arcTangent of the ratio of x and y.
  55.     }
  56. {}FUNCTION  D_to_R(x:REAL):REAL;
  57.     { Returns an angle converted to normalized radians (0 - 2π)
  58.     }
  59. {}FUNCTION  Exponent(x,n : REAL; VAR err:BYTE) : REAL;
  60.     { raise x to the n'th power (with error checking)
  61.     }
  62. {}FUNCTION  Gate(x,cntr,wide:REAL):REAL;
  63.     { Provides a rectangular pulse centered around CNTR of WIDE width.
  64.     }
  65. {}FUNCTION  Gaussian(x,cntr,variance:REAL):REAL;
  66.     { provides a normal pulse centered around CNTR (mean) with VARIANCE.
  67.     }
  68. {}FUNCTION  Log2(x:REAL):REAL;
  69.     { Returns Base 2 Log of x.
  70.     }
  71. {}FUNCTION  Log10(x:REAL):REAL;
  72.     { Returns Base 10 Log of x.
  73.     }
  74. {}PROCEDURE Normalize(value:REAL; VAR mantissa:REAL; VAR exponent:INTEGER);
  75.     { compute the (1<mantissa<10) and exponent of value
  76.     }
  77. {}FUNCTION  OrderOf(v:REAL):INTEGER;
  78.     { Returns the order of magnitude of v
  79.     }
  80. {}FUNCTION  PerCent_to_R(x:REAL):REAL;
  81.     { Converts percent slope to radians
  82.     }
  83. {}FUNCTION  Power( x, n : REAL ) : REAL;
  84.     { Returns X raised to the n'th power (no error checking)
  85.     }
  86. {}FUNCTION  R_to_D(x:REAL):REAL;
  87.     { Returns an angle converted to normalized degrees (0 - 360)
  88.     }
  89. {}FUNCTION  R_to_PerCent(x:REAL):REAL;
  90.     { Converts radians to percent slope
  91.     }
  92. {}FUNCTION  Sign( x:REAL ):REAL;
  93.     { Returns the sign of x, ie., -1, +1, or 0
  94.     }
  95. {}FUNCTION  Sinc(x,hite,freq:REAL):REAL;
  96.     { Sin(πƒx) / (πƒx)
  97.     }
  98. {}FUNCTION  Step(x:REAL):REAL;
  99.     { unit step function
  100.     }
  101. {}FUNCTION  Tan(x:REAL):REAL;
  102.     { Returns the Tangent of a value
  103.     }
  104. {}FUNCTION  Triangle(x,cntr,wide:REAL):REAL;
  105.     { provide a triangular pulse centered at CNTR of WIDE width.
  106.     }
  107.  
  108.      {====================================================================}
  109.  
  110. IMPLEMENTATION
  111. CONST
  112.   epsilon: REAL    = 1.0e-9;
  113.  
  114.  
  115. {}FUNCTION ArcCos(x:REAL):REAL;
  116. {+H
  117. ---------------------------------------------------------------------------
  118.   Purpose     - Returns the ArcCosine of x.
  119.  
  120.   Revised     - 1988.0202 (KSB) Vers:0.00; Wrote it.
  121.               - 1988.1026 (KSB) Vers:0.10; Fixed for small values and for ±1.-
  122. ---------------------------------------------------------------------------}
  123.   VAR
  124.     t    : REAL;
  125.   BEGIN
  126.     IF Abs(x) < epsilon THEN BEGIN        { if x nearly zero }
  127.       ArcCos := Pi/2.0;
  128.     END ELSE
  129.     IF (1.0 + x) < epsilon THEN
  130.          { if x = -1 }
  131.       ArcCos := Pi
  132.     ELSE
  133.     IF (x - 1.0) < epsilon THEN
  134.     { if x =  1 }
  135.       ArcCos := 0.0
  136.     ELSE BEGIN                          { if x is anything else }
  137.       t := ArcTan(Sqrt(1.0-Sqr(x))/x);
  138.       IF t < 0.0 THEN
  139.         ArcCos := t + Pi
  140.       ELSE
  141.         ArcCos := t;
  142.     END {IF};
  143. {}END {ArcCos};
  144.  
  145.  
  146.  
  147.  
  148. {}FUNCTION ArcSin(x:REAL):REAL;
  149. {+H
  150. ---------------------------------------------------------------------------
  151.   Purpose     - Returns the ArcSin of x.
  152.  
  153.   Revised     - 1988.0202 (KSB) Wrote it.
  154. ---------------------------------------------------------------------------}
  155.   BEGIN
  156.     IF (1.0 - Abs(x)) < 1.0e-9 THEN BEGIN
  157.       IF x < 0 THEN
  158.         ArcSin := -Pi/2.0
  159.       ELSE
  160.         ArcSin :=  Pi/2.0;
  161.     END ELSE
  162.       ArcSin := ArcTan(x/Sqrt(1.0-Sqr(x)));
  163. {}END {ArcSin};
  164.  
  165.  
  166.  
  167.  
  168. {}FUNCTION ArcTan2( y,x : REAL ) : REAL;
  169. {+H
  170. ---------------------------------------------------------------------------
  171.   Purpose     - Returns the ArcTangent of x.
  172.  
  173.   Revised     - 1988.0202 (KSB) Wrote it.
  174. ---------------------------------------------------------------------------}
  175.   CONST
  176.     eps  = 1.0e-4;
  177.   VAR
  178.     ix,iy: BYTE;
  179.   BEGIN
  180.     IF x < -eps THEN
  181.       ix := 255
  182.     ELSE
  183.     IF x > eps THEN
  184.       ix := 1
  185.     ELSE
  186.       ix := 0;
  187.  
  188.     IF y < -eps THEN
  189.       iy := 255
  190.     ELSE
  191.     IF y > eps THEN
  192.       iy := 1
  193.     ELSE
  194.       iy := 0;
  195.  
  196.     CASE ix OF
  197.       255 :
  198.       CASE iy OF
  199.         255 : ArcTan2 := Pi + ArcTan( Abs(y) / Abs(x) );
  200.         0 : ArcTan2 := Pi;
  201.         1 : ArcTan2 := Pi - ArcTan( Abs(y) / Abs(x) );
  202.       END {CASE};
  203.       0 :
  204.       CASE iy OF
  205.         255 : ArcTan2 := 3*Pi/2;
  206.         0 : ArcTan2 := 0;
  207.         1 : ArcTan2 := Pi/2;
  208.       END {CASE};
  209.       1 :
  210.       CASE iy OF
  211.         255 : ArcTan2 := 2*Pi - ArcTan( Abs(y) / Abs(x) );
  212.         0 : ArcTan2 := 0;
  213.         1 : ArcTan2 := ArcTan( Abs(y) / Abs(x) );
  214.       END {CASE};
  215.     END {CASE};
  216. {}END {ArcTan2};
  217.  
  218.  
  219.  
  220.  
  221. {}FUNCTION Log2(x:REAL):REAL;
  222. {+H
  223. ---------------------------------------------------------------------------
  224.   Purpose     - Returns the Base 2 Logrithm of x.
  225.  
  226.   Revised     - 1988.0202 (KSB) Wrote it.
  227. ---------------------------------------------------------------------------}
  228.   BEGIN
  229.     Log2 := Ln(x)/Ln(2);
  230. {}END {Log2};
  231.  
  232.  
  233.  
  234.  
  235. {}FUNCTION Log10(x:REAL):REAL;
  236. {+H
  237. ---------------------------------------------------------------------------
  238.   Purpose     - Returns the Base 10 Logrithm of x.
  239.  
  240.   Revised     - 1988.0202 (KSB) Wrote it.
  241.               - 1988.0929 (KSB) Added low value trapping
  242. ---------------------------------------------------------------------------}
  243.   BEGIN
  244.     IF x < 1.0e-12 THEN
  245.       x := 1.0e-12;
  246.     Log10 := Ln(x)/Ln(10.0);
  247. {}END {Log10};
  248.  
  249.  
  250.  
  251.  
  252. {}FUNCTION Tan(x:REAL):REAL;
  253. {+H
  254. ---------------------------------------------------------------------------
  255.   Purpose     - Returns the Tangent of x.
  256.  
  257.   Revised     - 1988.0202 (KSB) Wrote it.
  258. ---------------------------------------------------------------------------}
  259.   BEGIN
  260.     IF Abs(Cos(x)) < epsilon THEN
  261.       Tan := Pi/2
  262.     ELSE
  263.       Tan := Sin(x)/Cos(x);
  264. {}END {Tan};
  265.  
  266.  
  267.  
  268.  
  269. {}FUNCTION D_to_R( x : REAL ) : REAL;
  270. {+H
  271. ---------------------------------------------------------------------------
  272.   Purpose     - Returns an angle in radians.
  273.  
  274.   Revised     - 1988.0303 (KSB) Wrote it.
  275. ---------------------------------------------------------------------------}
  276.   BEGIN
  277.     D_to_R := 2.0 * Pi * Frac(x / 360.0);
  278. {}END {D_to_R};
  279.  
  280.  
  281.  
  282.  
  283. {}FUNCTION R_to_D( x : REAL ) : REAL;
  284. {+H
  285. ---------------------------------------------------------------------------
  286.   Purpose     - Returns an angle in degrees.
  287.  
  288.   Revised     - 1988.0303 (KSB) Wrote it.
  289. ---------------------------------------------------------------------------}
  290.   BEGIN
  291.     R_to_D :=  360.0 * Frac(x / (2.0*Pi));
  292. {}END {R_to_D};
  293.  
  294.  
  295.  
  296.  
  297. {}FUNCTION  PerCent_to_R(x:REAL):REAL;
  298. {+H
  299. ---------------------------------------------------------------------------
  300.   Purpose     - Returns radians for an angle expressed as % slope.
  301.  
  302.   Revised     - 1988.0901 (KSB) Wrote it.
  303. ---------------------------------------------------------------------------}
  304.   BEGIN
  305.     PerCent_to_R := ArcTan(x/100.0);
  306. {}END {PerCent_to_R};
  307.  
  308.  
  309.  
  310.  
  311. {}FUNCTION  R_to_PerCent(x:REAL):REAL;
  312. {+H
  313. ---------------------------------------------------------------------------
  314.   Purpose     - Returns % slope for an angle expressed as radians.
  315.  
  316.   Revised     - 1988.0901 (KSB) Wrote it.
  317. ---------------------------------------------------------------------------}
  318.   BEGIN
  319.     R_to_PerCent := Tan(x)*100.0;
  320. {}END {R_to_PerCent};
  321.  
  322.  
  323.  
  324.  
  325. {}FUNCTION Power( x, n : REAL ) : REAL;
  326. {+H
  327. ---------------------------------------------------------------------------
  328.   Purpose     - Raise X to the n'th power, Xⁿ (exponentiation).
  329. ---------------------------------------------------------------------------}
  330.   BEGIN
  331.     {Compute power using logarithms}
  332.     power := Exp( n * Ln( x ) );
  333. {}END {Power};
  334.  
  335.  
  336.  
  337.  
  338. {}FUNCTION OrderOf(v:REAL):INTEGER;
  339. {+H
  340. ---------------------------------------------------------------------------
  341.   Purpose     - Determines the order of magnitude of the argument
  342.  
  343.   Revised     - 1988.0708 (KSB) based on GetOrder.
  344. ---------------------------------------------------------------------------}
  345.   VAR
  346.     pwr  : INTEGER;
  347.   BEGIN
  348.     v := Abs(v);
  349.  
  350.     IF v = 0.0 THEN BEGIN
  351.       OrderOf := 0;
  352.       Exit;
  353.     END {IF};
  354.  
  355.     pwr := 0;
  356.     WHILE v >= 10.0 DO BEGIN
  357.       v := v/10.0;
  358.       Inc(pwr);
  359.     END {WHILE};
  360.  
  361.     WHILE (v < 1.0) DO BEGIN
  362.       v := v*10.0;
  363.       Dec(pwr);
  364.     END {WHILE};
  365.  
  366.     OrderOf := pwr;
  367. {}END {OrderOf};
  368.  
  369.  
  370.  
  371.  
  372. {}FUNCTION Sign( x:REAL ):REAL;
  373. {+H
  374. ---------------------------------------------------------------------------
  375.   Purpose     - Returns the sign of x, ie., -1, +1, or 0
  376. ---------------------------------------------------------------------------}
  377.   BEGIN
  378.     IF x < 0.0 THEN
  379.       Sign := -1.0
  380.     ELSE
  381.     IF x > 0.0 THEN
  382.       Sign := 1.0
  383.     ELSE
  384.       Sign := 0.0;
  385. {}END {Sign};
  386.  
  387.  
  388.  
  389.  
  390. {}FUNCTION Sinc(x,hite,freq:REAL):REAL;
  391. {+H
  392. ---------------------------------------------------------------------------
  393.   Purpose     - Returns sinc(x)
  394. ---------------------------------------------------------------------------}
  395.   BEGIN
  396.     IF Abs(x) < epsilon THEN
  397.       Sinc := hite
  398.     ELSE
  399.     IF Abs(freq) < epsilon THEN
  400.       Sinc := 0.0
  401.     ELSE
  402.       Sinc := hite * Sin(Pi*x*freq) / (Pi*x*freq);
  403. {}END {Sinc};
  404.  
  405.  
  406.  
  407.  
  408. {}FUNCTION Step(x:REAL):REAL;
  409. {+H
  410. ---------------------------------------------------------------------------
  411.   Purpose     - Provides unit step function
  412. ---------------------------------------------------------------------------}
  413.   BEGIN
  414.     IF x <= 0.0 THEN
  415.       Step := 0.0
  416.     ELSE
  417.       Step := 1.0;
  418. {}END {Step};
  419.  
  420.  
  421.  
  422.  
  423. {}FUNCTION Gate(x,cntr,wide:REAL):REAL;
  424. {+H
  425. ---------------------------------------------------------------------------
  426.   Purpose     - returns a rectangular pulse centered at CNTR of WIDE width.
  427. ---------------------------------------------------------------------------}
  428.   BEGIN
  429.     IF (x < (cntr-0.5*wide)) OR
  430.        (x > (cntr+0.5*wide)) THEN
  431.       Gate := 0.0
  432.     ELSE
  433.       Gate := 1.0;
  434. {}END {Gate};
  435.  
  436.  
  437.  
  438.  
  439. {}FUNCTION Triangle(x,cntr,wide:REAL):REAL;
  440. {+H
  441. ---------------------------------------------------------------------------
  442.   Purpose     - returns a triangular pulse centered at CNTR of WIDE width.
  443. ---------------------------------------------------------------------------}
  444.   BEGIN
  445.     IF (x < (cntr-wide)) OR
  446.        (x > (cntr+wide)) THEN
  447.       Triangle := 0.0
  448.     ELSE
  449.     IF Abs(wide) < epsilon THEN
  450.       Triangle := 1.0
  451.     ELSE
  452.       Triangle := (1.0 - Abs(x - cntr)/wide);
  453. {}END {Triangle};
  454.  
  455.  
  456.  
  457.  
  458. {}FUNCTION Gaussian(x,cntr,variance:REAL):REAL;
  459. {+H
  460. ---------------------------------------------------------------------------
  461.   Purpose     - returns a gaussian (normally distributed) pulse centered at
  462.                 CNTR (mean) with VARIANCE.
  463. ---------------------------------------------------------------------------}
  464.   BEGIN
  465.     IF variance > epsilon THEN
  466.       Gaussian := Exp(-Sqr(x-cntr)/(2*Sqr(variance))) /
  467.       Sqrt( 2*Pi*variance )
  468.     ELSE
  469.       Gaussian := 1.0;
  470. {}END {Gaussian};
  471.  
  472.  
  473.  
  474.  
  475. {}FUNCTION Exponent(x,n : REAL; VAR err:BYTE) : REAL;
  476. {+H
  477. ---------------------------------------------------------------------------
  478.   Purpose     - Raise X to the n'th power, Xⁿ (exponentiation).
  479. ---------------------------------------------------------------------------}
  480.   CONST
  481.     epsilon   = 0.000001;{tolerance to regard a number as an integer }
  482.   VAR
  483.     i    : INTEGER;
  484.     p    : REAL;
  485.   BEGIN
  486.     err := 0;
  487.     IF Abs(n - Round(n)) < epsilon THEN BEGIN   { x raised to the n }
  488.       p := 1.0;
  489.       IF n > 0.0 THEN
  490.         FOR i := 1 TO Round(n) DO
  491.           p := p*x
  492.           ELSE
  493.           IF x = 0.0 THEN
  494.             p := 0
  495.           ELSE
  496.             FOR i :=  - 1 DOWNTO Round(n) DO
  497.               p := p/x;
  498.  
  499.       Exponent := p;
  500.     END ELSE
  501.     IF x > 0.0 THEN BEGIN
  502.       p := n*Ln(x);
  503.       IF (p > 87) THEN
  504.         p := 87;
  505.       IF (p <  - 87) THEN
  506.         p :=  - 87;
  507.       Exponent := Exp(p);
  508.     END ELSE
  509.     IF Abs(x) < epsilon THEN
  510.       Exponent := 0
  511.     ELSE
  512.       err := 1;
  513. {}END {Exponent};
  514.  
  515.  
  516.  
  517.  
  518. {}FUNCTION Approx(a,b:REAL):REAL;
  519. {+H
  520. ---------------------------------------------------------------------------
  521.   Purpose     - A is an approximation value, 0<A<1. B is the value-
  522.                 to be adjusted. If the fractional part of B > A
  523.                 then round B up to the nearest whole number, else
  524.                 truncate B.
  525. ---------------------------------------------------------------------------}
  526.   BEGIN
  527.     IF (a>=1.0) OR (a<=0) THEN
  528.       Approx := Round(b)
  529.     ELSE BEGIN
  530.       IF Frac(b)>a THEN
  531.         Approx := Round(b)
  532.       ELSE
  533.         Approx := Trunc(b);
  534.     END {IF};
  535. {}END {Approx};
  536.  
  537.  
  538.  
  539.  
  540. {}PROCEDURE Normalize(value:REAL; VAR mantissa:REAL; VAR exponent:INTEGER);
  541. {+H
  542. ---------------------------------------------------------------------------
  543.   Purpose     - Given a real VALUE compute its (1 < MANTISSA < 10) & EXPONENT-
  544. ---------------------------------------------------------------------------}
  545.   BEGIN
  546.     exponent := 0;
  547.     mantissa := value;
  548.     REPEAT
  549.       IF (mantissa >= 10.0) THEN BEGIN
  550.         mantissa := mantissa / 10.0;
  551.         Inc(exponent);
  552.       END ELSE
  553.       IF (mantissa < 1.0) THEN BEGIN
  554.         mantissa := mantissa * 10.0;
  555.         Dec(exponent);
  556.       END {IF};
  557.     UNTIL (mantissa >= 1.0) AND (mantissa < 10.0);
  558. {}END {Normalize};
  559.  
  560.  
  561.  
  562.  
  563. END {UNIT}.
  564.